{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e352c272-9316-478d-8959-3cd643349329",
   "metadata": {},
   "source": [
    "Chapter 07\n",
    "\n",
    "# 初等行变换求解线性方程组\n",
    "《线性代数》 | 鸢尾花书：数学不难"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "447d4e7d-3744-42c8-b514-4f3d2795ee7d",
   "metadata": {},
   "source": [
    "该代码实现了高斯消元法（Gaussian Elimination）来求解线性方程组 $Ax = b$，其中 $A$ 是系数矩阵，$b$ 是右端项向量，$x$ 是待求解向量。高斯消元的核心思想是通过行变换将增广矩阵化为上三角矩阵，然后通过回代（Back Substitution）求解未知数。下面从数学角度详细描述代码执行的数学过程。\n",
    "\n",
    "---\n",
    "\n",
    "### **1. 形成增广矩阵**\n",
    "首先，将线性方程组 $Ax = b$ 组装成增广矩阵：\n",
    "$$\n",
    "[A | b] =\n",
    "\\begin{bmatrix}\n",
    "a_{11} & a_{12} & \\cdots & a_{1n} & | b_1 \\\\\n",
    "a_{21} & a_{22} & \\cdots & a_{2n} & | b_2 \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots & | \\vdots \\\\\n",
    "a_{n1} & a_{n2} & \\cdots & a_{nn} & | b_n\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "其中，$A$ 是 $n \\times n$ 的系数矩阵，$b$ 是 $n \\times 1$ 的列向量。代码通过 `np.hstack([A, b.reshape(-1, 1)])` 来构造这个增广矩阵。\n",
    "\n",
    "---\n",
    "\n",
    "### **2. 高斯消元**\n",
    "高斯消元的目标是通过一系列初等行变换，将增广矩阵转换为上三角矩阵：\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "1 & * & * & \\cdots & * & | * \\\\\n",
    "0 & 1 & * & \\cdots & * & | * \\\\\n",
    "0 & 0 & 1 & \\cdots & * & | * \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots \\\\\n",
    "0 & 0 & 0 & \\cdots & 1 & | *\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "在此过程中，主要涉及以下步骤：\n",
    "\n",
    "#### **(a) 选取主元并交换行**\n",
    "在第 $i$ 步时，代码使用 `np.argmax(np.abs(Ab[i:, i])) + i` 选择当前列中绝对值最大的元素作为主元，并交换当前行 $i$ 和主元所在行 `max_row`：\n",
    "$$\n",
    "\\text{若 } |a_{pi}| = \\max\\{|a_{ji}|, j \\geq i\\}, \\text{则交换 } i \\text{ 行和 } p \\text{ 行}\n",
    "$$\n",
    "这一操作可以减少数值误差，提高计算稳定性。\n",
    "\n",
    "#### **(b) 归一化主元**\n",
    "将当前主元归一化为 1，即对当前行的所有元素进行除法运算：\n",
    "$$\n",
    "\\frac{A[i, j]}{A[i, i]}, \\quad \\forall j \\geq i\n",
    "$$\n",
    "代码执行：\n",
    "```python\n",
    "Ab[i] = Ab[i] / Ab[i, i]\n",
    "```\n",
    "\n",
    "#### **(c) 消去主元下方元素**\n",
    "对所有 $j > i$ 的行执行消元：\n",
    "$$\n",
    "\\text{令 } m_{ji} = \\frac{A[j,i]}{A[i,i]}, \\quad A[j, k] = A[j, k] - m_{ji} A[i, k] \\quad \\forall k \\geq i\n",
    "$$\n",
    "等价于用第 $i$ 行的倍数消去第 $i$ 列下方的元素：\n",
    "```python\n",
    "Ab[j] = Ab[j] - Ab[j, i] * Ab[i]\n",
    "```\n",
    "最终，增广矩阵被化为上三角矩阵：\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "1 & a_{12}' & a_{13}' & \\cdots & a_{1n}' & | b_1' \\\\\n",
    "0 & 1 & a_{23}' & \\cdots & a_{2n}' & | b_2' \\\\\n",
    "0 & 0 & 1 & \\cdots & a_{3n}' & | b_3' \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots \\\\\n",
    "0 & 0 & 0 & \\cdots & 1 & | b_n'\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "---\n",
    "\n",
    "### **3. 回代求解**\n",
    "已得到上三角矩阵，接下来使用回代法计算 $x$：\n",
    "对于最后一行：\n",
    "$$\n",
    "x_n = b_n'\n",
    "$$\n",
    "对前面的行，逐步计算：\n",
    "$$\n",
    "x_i = b_i' - \\sum_{j=i+1}^{n} a_{ij}' x_j, \\quad \\forall i = n-1, n-2, \\dots, 1\n",
    "$$\n",
    "代码执行：\n",
    "```python\n",
    "x[i] = Ab[i, -1] - np.dot(Ab[i, i+1:n], x[i+1:n])\n",
    "```\n",
    "最终得到解向量：\n",
    "$$\n",
    "x = \\begin{bmatrix} x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "---\n",
    "\n",
    "### **4. 示例计算**\n",
    "设：\n",
    "$$\n",
    "A = \\begin{bmatrix} 1 & 3 & 6 \\\\ 2 & 7 & 14 \\\\ 0 & 2 & 5 \\end{bmatrix}, \\quad\n",
    "b = \\begin{bmatrix} 25 \\\\ 58 \\\\ 19 \\end{bmatrix}\n",
    "$$\n",
    "执行高斯消元后得到：\n",
    "$$\n",
    "\\begin{bmatrix} 1 & 3.5 & 7 & | 29 \\\\ 0 & 1 & 2.5 & | 9.5 \\\\ 0 & 0 & 1 & | 3 \\end{bmatrix}\n",
    "$$\n",
    "回代计算：\n",
    "$$\n",
    "x_3 = 3, \\quad x_2 = 9.5 - 2.5 \\times 3 = 2, \\quad x_1 = 29 - 3.5 \\times 2 - 7 \\times 3 = 1\n",
    "$$\n",
    "最终解为：\n",
    "$$\n",
    "x = \\begin{bmatrix} 1 \\\\ 2 \\\\ 3 \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "---\n",
    "\n",
    "### **5. 结论**\n",
    "代码实现了高斯消元法，先通过选取主元、行交换、归一化和消元操作，将增广矩阵转化为上三角形式，再使用回代求解线性方程组。该方法适用于求解具有唯一解的方程组，且数值稳定性较好。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "919eda14-8eb8-4e80-be0c-9dcde7dc6537",
   "metadata": {},
   "source": [
    "## 初始化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5d1a3ec3-302e-4963-8108-c24ccfc2a5de",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b518d9c3-5fc8-402d-ac88-82740a8fe79e",
   "metadata": {},
   "source": [
    "## 自定义函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d7bd2854-8517-46c8-9d24-a59c8ff6ac6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss_elimination(A, b):\n",
    "    # 确保使用浮点数计算\n",
    "    A = A.astype(float)  \n",
    "    b = b.astype(float)\n",
    "    n = len(b)\n",
    "\n",
    "    # 构造增广矩阵 [A|b]\n",
    "    Ab = np.hstack([A, b.reshape(-1, 1)])\n",
    "\n",
    "    # 原始矩阵A > 上三角\n",
    "    for i in range(n):\n",
    "        # 1) 行交换，选择绝对值最大为主元\n",
    "        # 主元：主对角线元素\n",
    "        max_row = np.argmax(np.abs(Ab[i:, i])) + i\n",
    "        Ab[[i, max_row]] = Ab[[max_row, i]]  \n",
    "        # 交换行, 行 i <> 行 max_row\n",
    "\n",
    "        # 2) 行数乘，归一化主元，主元变为 1\n",
    "        Ab[i] = Ab[i] / Ab[i, i]\n",
    "\n",
    "        # 3) 行倍加，消元，使主元下面所有行的 i列变为 0\n",
    "        for j in range(i+1, n):\n",
    "            # Ab[j] -= Ab[j, i] * Ab[i]\n",
    "            Ab[j] = Ab[j] - Ab[j, i] * Ab[i]\n",
    "\n",
    "    # 回代求解，上三角矩阵 > 单位矩阵\n",
    "    x = np.zeros(n)\n",
    "    for i in range(n-1, -1, -1):\n",
    "        x[i] = Ab[i, -1] - np.dot(Ab[i, i+1:n], x[i+1:n])\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32489668-e8b7-403d-9129-20d81944dcb7",
   "metadata": {},
   "source": [
    "## 定义矩阵 A 和向量 b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ee59f1e3-d975-4933-a96b-f33c6d6dbca4",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1, 3, 6], \n",
    "              [2, 7, 14],\n",
    "              [0, 2, 5]])\n",
    "b = np.array([25, 58, 19])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b64a3dbc-b6ff-429e-93bb-aef86045331f",
   "metadata": {},
   "source": [
    "## 求解 Ax = b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "98deb34b-9b04-4459-918d-e92e28a1cbf3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 2., 3.])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = gauss_elimination(A, b)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3b35558-e980-4a4d-ba1b-a4044700c3e0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "16650792-deab-4d0b-a223-2a6b7f41a5e1",
   "metadata": {},
   "source": [
    "作者\t**生姜DrGinger**  \n",
    "脚本\t**生姜DrGinger**  \n",
    "视频\t**崔崔CuiCui**  \n",
    "开源资源\t[**GitHub**](https://github.com/Visualize-ML)  \n",
    "平台\t[**油管**](https://www.youtube.com/@DrGinger_Jiang)\t\t\n",
    "\t\t[**iris小课堂**](https://space.bilibili.com/3546865719052873)\t\t\n",
    "\t\t[**生姜DrGinger**](https://space.bilibili.com/513194466)  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:base] *",
   "language": "python",
   "name": "conda-base-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
